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Abstract. We present an algorithm developed particularly to detect gravitationally lensed arcs in clusters of galaxies. This 
algorithm is suited for automated surveys as well as individual arc detections. New methods are used for image smoothing and 
source detection. The smoothing is performed by so-called anisotropic diffusion, which maintains the shape of the arcs and 
does not disperse them. The algorithm is much more efficient in detecting arcs than other source finding algorithms and the 
detection by eye. 
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1. Introduction 

Gravitational lensing has turned out to be a universal tool for 
very different astrophysical applications. In particular, lensing 
by massive clusters of galaxies is extremely useful for cos- 
mology. The measurement of various properties of the mag- 
nified and distorted images of background galaxies ("arcs and 
arclets") provides information on the cluster as well as on the 
population of faint and distant galaxies. Many of t hese distant 
backgr ound galaxies (up to redshifts of z ^ 5 , iFranx et all 
, (119971) ^ could not be studied with the largest telescopes if they 
were not magnified by the gravitational lensing effect. Some 
, of these distant gala xies are particularly useful for the s tudy 
' of galaxy evolution llSeitz et al Jll 9*9^ IPettini et a\\\20()(i . As 
these background galaxies are free of selection effects, be- 
cause they lie serendipitously behind massive clusters, they 
are ideal targets for a population study of distant galaxies 
llFortet alJll997l) . Gravitationally lensed arcs also provide a 
way to measure the total mass and the dark matter in clus - 
ters llFort & Mellielll994 IWambsgansslll998t iMemejl 19991) . 
As galaxy clusters can be considered to be fair samples of 
the universal mass fractions, such determinations probe cos- 
mological parameters like Qtot, ^^matter, and Qbaryon- A third, 
very important application of gravitational lensing in clus- 
ters is the determination of the frequency of arcs (=arc statis- 
tics). This is a strong criterion in order to distinguish be- 
tween different cosmological models (Bartelmann et al. 1998; 
iKaufmann & Strau mann 2000J. Therefore detections of gravi- 
tationally lensed arcs are of high importance for astrophysics 
and cosmology. 
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Ideal for cosmological studies are systematic searches. A 
successful arc search was perfor med with the X-ray lumi- 
nous cluster sample of the EMSS jLuppino et al.lfl99^ . More 
searches are under way, which not only cover larger areas than 
the previous survey, but they also go much deeper, i.e. fainter 
galaxies can be detected. 

The first ar c s we re detected only in 1986 
( iLvnds & PetrosianI Il986t 'Soucail e t al.l 1 19871) because 
they are very faint and very thin structures. Under non-ideal 
observing conditions (e.g. bad seeing) they are easily dispersed 
and disappear into the background. Even with ideal conditions 
they are not easy to detect because they are often just above 
the background level. In order to remove the noise and make 
faint structures better visible usually smoothing is applied. 
Unfortunately in the case of such thin structures as arcs the 
smoothing procedure often leads to a dispersion of the few 
photons so that the arcs are difficult to detect at all. To prevent 
this dispersing we suggest an algorithm that automatically 
smooths only along the arcs and not perpendicular to them, 
so-called "anisotropic diffusion". The subsequently applied 
source finding procedure extracts all the information from 
the sources necessary to distinguish arcs from other sources 
(i.e. edge-on spirals or highly elliptical galaxies). This new 
algorithm is much more efficient in finding gravitationally 
lensed arcs than existing source detection algorithms, because 
it is optimized just for this purpose. 

In Sect. 2 the algorithm is explained with its four different 
steps. In Sect. 3 examples of detected arcs are presented. Sect. 4 
outlines the differences to existing source finding software and 
the advantages for arc detection. In Sect. 5 we draw conclusions 
on the applicabifity and usefulness of the new algorithm. 
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2. The algorithm 

We propose a four level strategy for numerical detection of arcs 
in astronomical data sets consisting in the successive realiza- 
tion of 

1 . histogram modification 

2. anisotropic diffusion filtering 

3. object finding 

4. selection of arcs. 

The algorithm described in detail below has been implemented 
in the programming language C. 

The image data are given by a 2D matrix of intensity values 
I. For the sake of simplicity of presentation we assume that the 
intensity matrix is of dimension N xN. The set of indexes 

P :={(/,;•), 1...A?) 

are referred to as pixels and identify the position of the recorded 
intensity /(/, j). 




Fig. 1. Distribution of pixel values of an astronomical test im- 
age. A typical parameter setting for a and b (low and high cut) 
is marked in gray. A good choice for parameter a is the maxi- 
mum of the distribution and can be easily computed , whereas a 
general optimum for parameter b can not be prescribed. In the 
considered example a good choice is: a = and b - 1. 



2.1. Histogram modification 

Astronomical image data contain objects on a variety of bright- 
ness scales. Frequently stars and galaxies show up relatively 
bright; arcs however are small elongated objects of only 
marginally higher intensity than the surrounding background. 
In order to detect such arcs it is necessary to correct for the 
dominance of extremely bright objects. This is done by his- 
togram modification. Here, we use a nonlinear transformation 

{0 if X < fl 
sif^J if xe[a,b] (1) 
1 ifx>b 

which is applied to the pixel intensities I{i,j), giving a new 
intensity matrix /"(/, j) = S{I(i, j)). The transformation maps 
the pixel intensity distribution, which originally varied from 
to the maximum pixel value, onto the interval [0, 1] by the use 
of a bijective transformation s from the interval [0, 1] to [0, 1]. 

The interval [a, b] specifies the level of intensities where 
arcs are to be detected. The lower bound a is considered the in- 
tensity value of the background. By analyzing several different 
astronomical data sets we have learned that a and b have to be 
chosen relatively close to each other for optimal visualization 
of arcs (cf. Fig.Q. 

In the interval [a, b] we have to distinguish between noise 
and real sources such as stars, galaxies and arcs. To ease this 
separation process we apply nonlinear intensity transforma- 
tions such as s(x) - ^/x, or alternatively s{x) = x. 

Some astronomical data set may contain arcs in dififerent in- 
tensity ranges. In this case choosing the value for the parameter 
b of about the intensity of the brightest arcs is appropriate. 

2.2. Anisotropic Diffusion Filtering 

2.2.1 . Introduction to diffusion processes 

By interpolation the scaled image /"(/, j) can by interpreted as 
a function in and is now denoted by M°(x,y), ix,y) e M^. 



By applying Gaussian convolution with a kernel K,(x,y) :- 
^ exp (^—^jf-^ depending on t one gets a smoothed image 

u(t,x,y) :- u^(x,y) * Ki(x,y) 

:- / u'^{x - x,y — y) K,(x,y) dxdy. 

The parameter t controls the amount of smoothing. A larger 
value of t corresponds to a higher filtering. 

In the following dru{t, x, y) denotes the first derivative of u 
with respect to f, dxxu{t, x, y) and dyyu{t, x, y) the second deriva- 
tives with respect to x and y and tS.u{t,x,y) :- dxxu(t,x,y) + 
dyyu(t, X, y) the Laplacian. 

It is well known that u{t,x,y) solves the diffusion (resp. 
heat) equation 

d,u(t, X, y) - Auit, x,y)^Q for all (x, y) e (2) 

with the initial condition 

m(0, X, y) = m"(jic, y) 

To simplify the notation we write u(t, x), x e instead of 
u(t, x,y). 

Eq. (|3 can be restricted to a rectangular domain Q c M.^, 
the domain of interest, typically the set of pixels where inten- 
sity information has been recorded. 

In order to achieve existence and uniqueness of the solution 
u(t, x), X € O it is necessary to prescribe boundary conditions 
such as 

dn 

— (f, x) = for all t e (0, oo) and x e 5Q, 
on 

where c?0 denotes the boundary of the domain Q, n is the outer 
normal vector on c?Q and ^ denotes the derivative in direction 
of n. 

Applying the diffusion process up to a fixed time T > 
smooths the given data and spurious noise is filtered. The 



F. Lenzen and S. Schindler and O. Scherzer: Automatic detection of arcs and arclets 



3 



parameter T defines the strength of the filtering process. Thus 
in the following we refer to T as the filter parameter. 

The disadvantage of Gaussian convolution is that edges in 
the filtered image u(T,x) are blurred and the allocation and 
detection of object borders is difficult. To settle this problem 
several advanced diffusion models have be en proposed in the 
literature (Weicke rIl998llCatte et alJl992h . 

In the next section we define the general model of a diffu- 
sion process and in Sect. 12. 2.31 we describe the specific model 
used in our algorithm. 

2.2.2. The general diffusion equation 

Anisotropic diffusion filtering consists in solving the time de- 
pendent differential equation 

d,u(t, x) - div (jD(x, u, Vu)'Vu(t, x)j = (3) 

up to a certain time T > 0. 

Here D(x, u, Vm) is the 2x2 diffusion matrix depending on 
X and u. 

We prescribe the same initial and similar boundary condi- 
tions as mentioned above. Setting D(x, u, Vm) = 1 results in the 
heat equation Q. 

Two classes of anisotropic diffusion models are considered 
in the literature: if D{x, u, Vm) is independent of u and Vm then 
Eq. Q is called linear anisotropic diffusion filtering, otherwise 
it is called nonlinear filtering. 

For a survey on the topic of diffusion filtering we refer to 

IWeickerttiggS.^ . 

2.2.3. Tine specific model 

In anisotropic models the diffusion matrix D is constructed to 
reflect the estimated edge structure. That is to prefer smoothing 
in directions along edges, or in other words edges are preserved 
or even enhanced and simultaneously spurious noise is filtered. 

Consequently this kind of filtering eases a subsequent edge- 
based object detection. 

In the following D will only depend on the gradient Vm, 
which reflects the edge structure of the image. 

To accomplish the diffusion matrix D(Vm) we note that 

Vm Vm""- 1 / du du\ 

vi := and vo '■- '■- — , 

|Vm| |Vm| |Vm| \dy dx J 

denote the directions perpendicular and parallel, respectively, 

to edges. (By v"" = | ) = | ) we denote the vector 

\V2J \-viJ 

perpendicular to v). 
By selecting 

D(Vu)^(vuv,)^^^^l''^^'^yvuV,f. 

with 

p(|Vm|) :- ^ ~ , with parameter K > 



the diffusion filtering method Q prefers filtering parallel to 
edges. 

Fig. 13 highlights the diffusion directions: arrows indicate 
the main directions of diffusion (vi, V2) and their thickness re- 
lates to the diffusion coefficient determining the strength of 
diffusion. Parallel to the edge the diffusion coefficient is con- 
stantly 1 (strong diffusion), where as the diffusion coefficient 
g(|VM|) in normal direction decreases rapidly (c.f. Figure|2ji as 
|Vm| increases (weak diffusion on edges). 

The dependence of g(|VM|) on |Vm| is controlled by pa- 
rameter K (c.f. Fig. O. We therefore refer to K as the edge 
sensitivity parameter. 

We use the following common variation of the diffusion 
matrix D: 

To limit the effect of noise the diffusion tensor is chosen 
to be dependent on the pre-smoothed image Mo-(f, x) - u{t, x) * 
Kcr obtained by Gaussian convolution with pre-filter parameter 
cr > 0. In the following to determine the diffusion matrix we 
exploit the gradient of of the pre-filtered image Ua-{T, .) instead 
of u(T, .). 

Let v'l', be the eigenvalues of the filtered structure tensor 
J^(x):= {K^*[(Wu^)(Wu„f])(x) . 

We refer to yu > as the pre-filter parameter for the structure 
tensor. 

Note that in the case ji - the eigenvectors of Jq are 
and 1^^. If cr and p are small, the effect of Gaussian filtering 
is negligible, and consequently v^^, \^ and vi, V2 refer to similar 
edge structures. 

The purpose of filtering of the structure tensor is to average 
information about edges in the image. 

Taking into account the approximation we are led to the 
following diffusion matrix 

D„.(Vm.) = «,v^)(^^I''o"'^I'^ J)(v^,v^)^ 

Besides the fact that D^ uCVug-) is less noise sensitive there 
are several advantages of using the approximation D^_o-(Vmo-) 
instead of D{yu): 

1. The numerical calculation of DiVu) is unstable; if p and cr 
are chosen appropriately D^^a-i^Ua-) can be determined in a 
stable way. 

2. For our particular application, in a neighborhood of arcs the 
vector field is nearly parallel to the arcs. Thus filtering 
is performed in arc orientation enhancing the geometrical 
structure of arcs. A side effect which is very useful for our 
particular application is that small gaps between elongated 
structures are closed, merging nearby objects. 

This is the anisotropic diffusion filtering method used in our 
numerical calculations. 

In order to solve the differential equation Q it is dis- 
cretized using s. finite element method in space and an implicit 
Euler method in time. (Within the finite element method the 
width of the quadratic elements is set to 1.) For each time 
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Fig. 2. Graph plot of function g(x) used for weakening the dif- 
fusion orthogonal to edges to archive edge enhancing. For pa- 
rameter K values = 0.1 resp. K - 0.001 are used. 



step the resulting system of linear equations is solved by a 
conjugate gradient method. For a survey on solving parabolic 
different i al equ ations with finite element methods we refer to 
Thomee ( 1984). The conjugate gradient method is discussed in 
Hanke-Bourgeois ( 1995. 20021) . 

Anisotropic diffusion filtering requires to select parameters 
T, K, cr and yU. 

- Filter parameter T defines the amount of filtering. 
Consequently for small values of T, spurious noise is still 
recognizable. If T is large, spurious noise is eliminated but 
also small details are blurred. The parameter T has to be 
adapted to each data set to be analyzed. 

- Edge sensitivity parameter K determines the dependency 
of the strength of diffusion orthogonal to edges on |Vmo-|. 
The smaller K is, the more edges are enhanced. Thus, if 
K is chosen too low, the diffusion may generate artifacts 
out from the noise. Using a too large value for K leads to a 
filtering similar to Gaussian convolution, i.e. the image gets 
blurred. 

The parameter can be calibrated using test data and remains 
unchanged for similar data. 

- Pre-filter parameter cr can be selected from calibrated data. 

- The pre-filter parameter for the structure tensor ji controls 
the local average information. It has to be adopted to the 
size of the object to be recovered. 



2.3. Object finding 

In this section we discuss a partitioning algorithm to sepa- 
rate diff'erent image data, i.e. disjoint subsets of connected ob- 
jects and background {partitions). The algorithm uses only the 
anisotropic diffusion filtered data «(■) := u{T, ■) and not the ini- 
tial data m". 

In order to save computational eff'ort we restrict our atten- 
tion to segment objects of interest, i.e. isolated objects exceed- 
ing a certain brightness. We search for local intensity maxima 




Fig. 3. Diffusion (smoothing) near edges: The thickness of the 
arrows indicate the strength of diffusion. Parallel to the edge a 
strong diffusion occurs whereas in orthogonal direction a weak 
diff'usion leads to an enhancing of the edge. Due to the averag- 
ing of the structure tensor these directions are also determine 
the diffusion in a surrounding area of the edge and in partic- 
ular at the vertices yielding a diffusion mainly parallel to the 
direction of elongation. 



exceeding a certain intensity c^in (referred to as the intensity 
threshold for detection). 

Each maximum serves as seed for the partitioning algo- 
rithm: Starting from the seed pixel the region to which this 
pixel belongs has to be determined. 

To outline this concept we use the following notation. For 
a given pixel p - {i, j) e f we denote by A^(/:i) the set of the 
eight neighboring pixels. 
The neighborhood of a set c P is the set 

N{R) ^{qe Nip) : peR] . 

For two pixels p,q e P a path connecting p to qis any sequence 
(P = Po,PuP2, ...pn^q) satisfying pi e N(pi-i), /=!...«. 

The partitioning algorithm consists of two loops. 

1. Set j-0 and Pq := P for initialization. 

2 . Update j ^ ; + 1 ■ 

A pixel p e Pj_\ is selected where the 
intensity of the filtered image u{p) attains 
a local maximum exceeding the intensity 
threshold for detection: u{p) > Cmin > 0. 
A numerical procedure for detection of local 
maxima is described in Sect. 12 . 3 . 2L 
If no such pixel can be found the 
algorithm is terminated and Pj-\ is called 
' 'background' ' . 

3. The second step is a region growing 
algori thm. 

(a) It is initialized with / - 0, the seed 
po ■- p and the region Rq :- {po}- For po a 
threshold parameter cthies > is selected, 
which is used to terminate the region 
growing algorithm. The determination 
of this parameter is explained in 
Sect. 12.3.11 . 
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(b) Update 

Rj consists of and the neighboring 
pixels N(Ri-i) n Pj^i satisfying that 

- the according intensity exceeds Cthres 
and 

- the intensity is smaller than the 
intensities of neighboring pixels in 

Ri-i . 

(c) The iteration is terminated when = 
and we denote R(po) := R(i). This region 
corresponds to an object. 

4. Updating Pj = Pj-i\R(po) and repeating steps [3] 
and [B], completely determines the algorithm. 
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Fig. 4. Segmentation of an elliptical object (the theoretical ob- 
ject boundary is indicated by the black ellipse, gray squares 
indicate pixels with high intensities / white squares indicate 
pixels with low intensities): Assuming that the region growing 
starts with the pixel numbered by zero the neighboring pixels 
with numbers 1 to 7 are successively added to the region un- 
til the intensities of the next pixels to be added (here: white 
colored ) fall below a certain threshold intensity. 



Fig. 13 illustrates the process of region growing for an 
elliptical objects. 

A detailed overview ov er segmentatio n methods is given in 
lRosenfeld&Kai3 ill 9931) or lSoillelll2003h . 

In the following we describe numerical procedures for cal- 
culation of local maxima and the determination of Cthres in the 
object finding algorithm. 

2.3.1. Determination ofcthres 

Let po be a local maximum within an arc, which is used as an 
initialization for the region growing algorithm. The anisotropic 
diffusion filtering method calculates the two eigenvalues and 
of the structure tensor 7^. Thinking of an arc as an elliptically 
shaped region and approximate the axes of the ellipse, 
v^, denote the principal, respectively cross-sectional axis. 
Let M_„, M-„+i, . . .uq, . . .u„ be the intensities along the cross- 
sectional axes. Fig.|5lshows such a typical intensity distribution 
along a cross-sectional axes of an object. The points on the 
cross sectional axis with maximal gradient are natural candi- 



u 




j- Po r 

Fig. 5. Typical intensity distribution along a cross-section 
through an object 

dates for object boundaries. In the discrete setting these points 
are 

f = arg max {|m; - m,+i|) , 

i=-l...-n 

f = argmax{|M/ - m,_i|) . 

As threshold intensity for the region growing algorithm we 
choose: 

Cthres ■— max{Mj- , Mj+ ) 

2.3.2. Detection of local maxima 

For the detection of local maxima which exceed the inten- 
sity threshold for detection, c^in, we proceed in a simil ar way 
as in the implementation of watershed algorithms (c.f. 'S toevI 
( 2000)). The strategy for finding a local maximum is to choose 
an initial pixel p and to look for a pixel q neighboring to p with 
higher intensity and then with p ^ q reinitialized to proceed 
iteratively until a local maximum is reached. 

Applying this procedure we have to deal with the case, that 
a local maximum may not be a single pixel but a connected set 
of pixels with the same intensity (a plateau). In the sequel, for 
simplicity of presentation, we also refer to a single pixel as a 
plateau. 

Taking this into account we use the following procedure: 

We use markers (+), (-), and (0) to denote if a pixel is a 
local maximum, not a local maximum, or if it is not yet con- 
sidered in the algorithm, respectively. Initially every pixel is 
marked by (0). 

1. Search for a pixel p marked with (8). If 
no such pixel exists, the algorithm is 
terminated. 

If the intensity u(p) lies below c^in mark p 
with (-) and repeat this step with another 
pixel . 

2 . Start segmentation of the plateau including 
pixel p by applying a region-growing 
algorithm. During this process we monitor 
the following occurrences: 
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(a) A pixel q neighboring the plateau is found 
with u{q) > u{p): Mark the pixels of the 
plateau with (-) , set p - q and repeat 
step 2 . 

(b) A pixel of the plateau is found which is 
already marked with (-) : Mark the pixels 
of the plateau with (-) and go to step 1. 

(c) All neighboring pixels of the plateau 
have less intensity: Mark the pixels of 
the plateau with (+) and go to step 1 . 

After finding a local maximal plateau we choose one pixel 
of the plateau as seed for the later object detection. 

Note that the identification of local maximal plateaus can 
not be carried out by investigating the first and second deriva- 
tives of the intensity function. 

Threshold c^in correlates to the parameters a and b. The 
latter should be chosen such that a is about the background 
intensity and b is about the intensity of the arcs to be detected. 
As the arc intensities are only very little higher than the noise 
amplitude, Cmin serves as a threshold between the noise level 
and arcs intensity range and is of high importance for the object 
detection. 

Choosing Cmin ~ guarantees detection of most objects 
but many of them may result from noise amplification. On the 
other hand a choice Cmin ~ 1 real arcs with intensity below Cmin 
are obviously not detected. 



2.4. Selection of Arc Candidates 
2.4.1 . Selection of Source Candidates 

Due to blurring effects on the CCD, in a neighborhood of bright 
stars and galaxies, the background shows up bright, amplify- 
ing the intensity of the noise. In these regions several local in- 
tensity maxima above the threshold c^in occur and regions are 
detected, which clearly are artificial and belong to the back- 
ground. Such regions are singled out by a comparison of the 
mean intensities along a cross section within an object and the 
background. To this end let 



^back ■ — 



Ui+ ^ Ui 

.i=2j-...(j--i) ;=a++i)...2j+ 



The artificial objects (e.g. resulting from the diffusion filter- 
ing of noise) satisfy the condition that Mobj - Wback is small. 
Consequently we further restrict our attention to objects sat- 
isfying 

Mobj - "back > Cniin2 > 

The parameter c^ina denotes the minimal difference of the ob- 
jects intensity from the background. We refer to it as the pa- 
rameter of minimal intensity difference from background. 



By taking into account only pixels with indices between 
2 and 2 7^, Wback averages the intensities in a small neighbor- 
hood of the object. 

We assume that the object is isolated, i.e. that there are no 
other objects in this neighborhood and Mback reflects the local 
average background intensity. 

2.4.2. Selection of Arcs 

The last step of our algorithm consists in selecting arcs from 
the detected objects. For deriving the following features of 
the objects the image resulting from applying the histogram 
modification is used. 

Let R be an object, consisting of the set of pixels 

{pi^(xlx^):ieN]QP 

The light intensity of such an object can be interpreted as 
mass resp. density and thus regarding the object as an rigid 
body we define the total intensity 

m ^ u^\pi) . 

ieN 



The object's center is given by 
1 

mc^ 



mc 

mc = I 2 



m ^-^ \ X; J 



We define the 2nd moment M as the 2x2 matrix with compo- 
nents 



The eigenvalues Ai > /I2 of M indicate the length and the thick- 
ness of the object. 

The eccentricity of R 

Ai - A2 A2 

ecc :— = 1 

^1 ^1 



is a measure of elongation of R (I.T ahnel2002h . We identify arcs 
as thin and elongated objects, which in mathematical terms re- 
quires that 

ecc > Cecc e [0, 1] and A2 < Ctwck ■ 

for given thresholds Cecc (eccentricity threshold and CtMck 
(thickness threshold). 

The detection and selection of arcs is controlled by the three 
parameters c^ina, c^cc and CtMck: 

- Cmin2 (minimal intensity difference from background) is 
used to select arcs influenced by bright structures such as 
galaxies and dominating stars, thus it becomes active only 
locally. Choosing Cinin2 = is in general adequate. 

- Cecc 6 [0, 1] and Cthkk are statistical parameters of the shape 
of an object, they are used to select elongated and thin struc- 
tures. Choosing Cecc = allows the selection of spheres 
with a maximal diameter related to cthick- 
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It has to be taken into account, that an object may contain 
several local maxima and the region growing procedure detects 
adjacent parts of this object. These parts can be merged into 
one object after the selection process easily . Note that our 
sequence of selection and merging prevents objects of different 
shape from being merged. 

In most cases the criteria described above are sufficient to 
detect arcs. The selection process can be refined by incorpo- 
rating a priori information, like for instance information on 
the center of mass of a gravitational lens (galaxy cluster). For 
a spherically symmetric ideal gravitational lens with center 
gc - {gc\,gC2), the arcs occur tangentially around the center. 
Let mc denote the center of mass of an arc. Then the vectors 
p - mc - gc (position relative to center) and (approx. direc- 
tion of elongation) have to be orthogonal, giving an additional 
selection criterion for arcs. However, note that user interaction 
is required to incorporate a priori information on the position 
of the gravitational lens. 

The algorithm might be also used for detection of strings, 
if additional post-processing steps are applied using alternative 
selection criteria which take into account alignment informa- 
tion (instead of shape information as above). 



3. Results 

3.1. Quality of detections 

The quaUty of the results provided by our algorithm depends 
mostly on noise variance present in the data. 

Arcs may not be detected (false negative detection) by the 
algorithm if their intensity is within the scale of the noise. As 
we have seen in Sect. 12. li the intensity range of the arcs is in 
general close to the background intensity. The choice of pa- 
rameter Cmin determines whether a weak structure is interpreted 
as background noise (ignored) or is segmented (feasible local 
maxima). Corresponding the choice of a high value for c,nin in- 
creases the risk of a false negative detection. 

Another case leading to a false negative detection is the 
joining of close-by objects. 

A false positive detection may occur if the noise level 
exceeds the intensity Cmin- In such case the edge enhancing 
anisotropic diffusion may recover structures present in the 
noise, which may be segmented as elongated objects after- 
wards. Thus the choice of a low value for c^in increases the 
risk of a false positive detection. 

In the neighborhood of bright sources the background 
intensity increases and the noise may lead to several local 
maxima above c^nin- As a result the risk of a false positive de- 
tection increases in these areas. To reduce these kinds of false 
positive detections we utilize the parameter Cmin2 prescribing 
the minimal intensity difference of a detection in comparison 
with its surrounding. Again the quality of detection depends 
on the adaption of Cinin2 to the noise variance. 



3.2. Test images 

To highlight the properties of our algorithm we applied the al- 
gorithm to three astronomical test images. The first data set 
is an image of size 2285 x 2388 pixels with intensity range 
[-8.49,700.49], the second and third are of size 2048 x 2048 
pixels with intensity ranges [0, 19559.8] and [0,9314.26], re- 
spectively. We plotted in Fig. ^ the histogram of the first data 
set, the histograms of the second and third test image look sim- 
ilar. 

Figs- m Eland |8l show galaxy clusters with gravitationally 
lensed arcs, which have been treated with histogram transfor- 
mation; we have used s(x) = y/x and (a, b) - (0, 1), (149, 200) 
and (141, 170), respectively. The histogram modification is al- 
ready useful to visualize arcs. 




Fig. 6. Detail of a VLT observation of the galaxy cluster 
RXJ1347-1145 (from Miralles et al., in prep.). The image has 
a size of 2285 x 2388 pixels. This is the first test image. 



3.3. Computational effort 

Fig. |9] shows the computation time' for anisotropic diffusion 
filtering, object finding (including the search for local maxima, 
segmentation and selection) and the total computational time 
in dependence of the size of the image data (number of pixels) 
in typical astronomical data sets. 

Numerically, the pre-filtering step is most expensive as a 
large system of linear equations has to be solved. The num- 
ber of iterations performed by the conjugate gradient solver in- 
creases slowly with growing data size. Thus the computational 

' Computed on an Intel Pentium IV with 1.5 GHz 
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Fig. 7. Detail of an HST observation of the galaxy cluster 
A1689. The image has a size of 2048 x 2048 pixels. This is 
the second test image. 




Fig. 8. Detail of an HST observation of the galaxy cluster 
A1689. The image has a size of 2048 x 2048 pixels. This is 
the third test image. 



effort is approximately linearly correlated to the number of pix- 
els (cf.Fig.|9}. 

During the detection of local maxima each pixel is invoked 
a fix number of times: 

1 . to check if it is already marked, 

2. to check if it has to be included into a plateau during the re- 
gion growing process (once for each of its eight neighbors), 

3. to mark it as being (or not being) a local maximum. 

Therefore the computational effort for detecting the local max- 
ima is linearly depending on the number of pixels. As the seg- 
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Fig. 9. Computation time versus image size for the different 
steps of the algorithm. 



mented area is small in comparison with the image size, the 
effort for object finding strongly depends on the number of ob- 
jects resp. local maxima. The same holds for the selection pro- 
cess. 

Overall the total computation effort grows linearly with the 
image size (cf. Fig.|9j. 

3.4. Anisotropic Diffusion Filtering 

Fig. [TO] shows the result of the anisotropic diffusion filtering 
for the first test data set. To examine the effect of enhancing 
elongated structures we zoom in the filtered image. Fig. ^2 
shows the histogram modified image, i.e. the cuts of the 
image are set to the steep region of the distribution of the 
pixel values (top), the Gaussian filtered image (middle) and 
the image filtered with anisotropic diffusion (bottom). Both 
filters were applied separately to the histogram modified image. 

The filter parameters were chosen to remove noise up to 
nearly the same signal-to-noise ratio ^ , which is about 6.1 % 
in the Gaussian filtered image and about 6.3 % in the image 
filtered with anisotropic diffusion. 

In comparison to Gaussian convolution, anisotropic fil- 
tering is able to preserve accurately the edges of the objects 
both of high and of weak intensities. Anisotropic diffusion 
therefore is ideal as preprocessing for object finding based on 
edge detection. 



3.5. Detected objects 

In the first test image (Fig. |6} after applying the histogram 
modification about four arcs can be recognized at first glance. 
These arcs are grouped around a center of the cluster of 



defined as SNR= 



Jn("°(v)) ''-V 



with «" the original data set and 



JpX^fi{x}-u(T..y)) dx 

ii(T, .) the filtered image using Gaussian convolution resp. anisotropic 
diffusion. 
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Fig. 10. Image of the cluster RXJ1347-1 145 with anisotropic 
diffusion filtering applied. Compared to the unsmoothed image 
in Fig. |6lthe noise is considerably reduced. Parameter setting 
for the filtering: Filter parameter T - 15, edge sensitivity K = 
0.0001, pre-filter parameter cr = 2 and pre-filter parameter for 
structure tensor /i = 9 



galaxies, which appears in the middle of the image. 

In Fig. 1121 shows selected objects with thickness A2 < 16 
and eccentricity ecc > 0.7. The eccentricity is color coded: 
green, yellow, and red corresponds to an eccentricity in the 
ranges [0.7,0.8], [0.8,0.9] and [0.9, 1]. The higher the eccen- 
tricity and the smaller the thickness are, the higher is the prob- 
ability that the detected object is an arc. Incorporating a priori 
information on the center of the gravitational lens, unreason- 
able arcs candidates can be filtered out further (Fig.ll3>. Table 
[Olists the coordinates of the detected objects with an eccentric- 
ity greater than or equal to 0.84 (referring to Fig. ll2t . Besides 
the arcs already mentioned the algorithm detects a significant 
amount of arc candidates which are not obviously recognizable 
to the naked eye. 

Figs. ll5landll6lshow the results of our algorithm applied 
to the second and third test image. 

4. Comparison with software for source extraction 

In this section we compare our algorithm with the soft- 
ware package "SExtra ctor" by E. Bertin ( Berlin 1994| 
Berlin & Arnoutslll996l) . SExtractor is a general purpose 
astronomical software for the extraction of sources such as 
stars and galaxies, while our software is particularly designed 
to extract gravitationally lensed arcs. Although the main 
areas of applications are rather different, several levels of 
implementation are similar, although quite different in details: 




Fig. 11. Zoom of Fig. |6] Top: histogram modified data, mid- 
dle: Gaussian filtered image (Kernel K^r with cr - 5), bottom: 
image filtered with anisotropic diffusion with the parameters: 
filter parameter T =15, edge sensitivity K = 0.0001, pre-fiJter 
parameter cr - 2 and pre-filter parameter for structure tensor 
jj. = 9. In the Gaussian filtered image (middle) the edges are 
not preserved well, i.e. the arcs get dispersed, while anisotropic 
diffusion (bottom) maintains the edges and reduces the noise at 
the same time. 

Table 1. List of objects detected in the first data set (cf. Figs.|6l 
and ll2> with eccentricity larger than 0.84, i.e. good arc candi- 
dates. Objects with an eccentricity larger than 0.9, i.e. objects 
1 to 6 are very good candidates. 





Mass center 


Size (Pixel) 


Eccentricity 


1 


(593,749) 


768 


0.978 


2 


(188,177) 


603 


0.956 


3 


(595,899) 


143 


0.943 


4 


(794,440) 


241 


0.942 


5 


(320,144) 


196 


0.928 


6 


(842,901) 


181 


0.901 


7 


(141,328) 


27 


0.876 


8 


(130,819) 


95 


0.871 


9 


(165, 61) 


52 


0.866 


10 


(223,871) 


86 


0.864 


11 


(488,103) 


65 


0.854 


12 


(151,659) 


32 


0.845 


13 


(677,260) 


44 


0.840 



SExtractor uses background estimation which in our pro- 
gram is performed by histogram modification. For Detection, 
SExtractor used Gaussian convolution filtering; this step 
relates to the object finding process. Deblending and filtering 
of deblending are related to our merging strategy described at 
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Fig. 12. Result of segmentation and selection applied to the 
first test image - cluster RXJ1347-1 145. Only the objects 
which are detected as being arcs with high probability are 
shown. Settings: intensity threshold for detection c,nin = 0.1, 
minimal intensity difference from background c„iin2 =0.1, ec- 
centricity threshold Cecc - 0.7 and thickness threshold Cthkk = 
16. 

Table 2. Parameter settings used for the three test data sets. 



Param. 


Meaning 


Image 1 


Image 2 


Image 


a 


low cut 





149 


141 


b 


high cut 


1 


200 


170 




intensity threshold 


0.1 


0.1 


0.3 




for detection 








T 


filter parameter 


15 


30 


20 


K 


edge sensitivity 


0.0001 


0.0001 


0.0001 


cr 


pre-filter parameter 


2 


2 


2 




pre-filter parameter 


9 


30 


15 




for structure tensor 








^min2 


minimal intensity 


0.1 


0.2 


0.2 




diff. from backgr. 








Cecc 


eccentricity threshold 


0.7 


0.7 


0.7 


Cthick 


thickness threshold 


16 


9 


9 



tion several other objects with only small elongations show up 
when applying SExtractor. To distinguish close-by objects 
we use two different colors. 

The first remarkable difference is that SExtractor in order 
to measure the objects total magnitudes uses a far lower seg- 
mentation threshold and thus it segments larger parts of bright 
objects than our algorithm as can be realized in columns one to 
three in Fig.^^ The evaluation of the objects elongation then 
depends on the segmented shape. 

The results of our algorithm show a more regular border due 
to the edge parallel diffusion and the edge enhancement. Using 
SExtractor one may choose a higher threshold for detection 
(DETECT_THRES) and yield a better shape of the segmented ob- 
jects. However, since this is a global threshold SExtractor 
tends to loose fainter objects. Our algorithm overcomes this 
problem by using a local adaptive threshold. 

The forth arc in Fig. 1141 segmented by SExtractor is an 
example for a failure of the deblending procedure. The arc is 
not separated from the nearby galaxy and since the composition 
of both sources is not much elongated it would be refused by 
the selection process. On the other hand tuning the parameter 
for deblending, which is supplied by SExtractor for this spe- 
cific arc allows a separate detection of both sources but leads 
to an undesired splitting of other objects. 

In the fifth arc our algorithm reveals a far larger part 
of the weak structure than SExtractor due to the use of 
anisotropic diffusion and detects its direction of elongation cor- 
rectly. SExtractor does also find an elongated part of this arc 
but the direction of elongation of this part differs from the ex- 
act direction. Thus our algorithm provides a better quality of 
detection. 

To summarize the new algorithm for detecting arcs pre- 
sented here has two main advantages: 

The method of filtering, e.g. anisotropic diffusion is chosen 
with regard to the kind of objects to be detected. The filtering 
process provides a closure of gaps in between objects as well 
as edge enhancing. 

The use of object dependent thresholds based on edge 
detection leads to an improved segmentation even of weak 
sources. A deblending procedure is not required. Moreover the 
new algorithm is able to detect close-by objects of different 
scale, for which the SExtractor ' s deblending procedure 
fails. 



the second last paragraph of the last section. 

To compare both algorithms we single out five specific arcs 
in the first test image. Fig. ll4lshows these objects as they are 
detected by the proposed algorithm (upper row) and as they are 
segmented by SExtractor (lower row). 

Concerning the SExtractor ' s results a possible adaption 
of the SExtractor ' s algorithm for detection of arcs would 
be to perform a selection process afterwards as described in 
Sect. l2.4.2l However we did not apply such a selection process 
to the SExtractor ' s output. Beside the arcs under investiga- 



5. Summary and Conclusions 

We proposed a new algorithm for the detection of gravitation- 
ally lensed arcs on CCD images of galaxy clusters, performing 
an edge-based object detection on the filtered image together 
with an automatic selection of arcs. 
The algorithm consists of several steps: 

1 . Histogram modification: We focus on the typical range of 
the arcs intensities. 

2. Filtering: we use anisotropic diffusion, which enhances thin 
and elongated objects such as arcs. 
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Fig. 13. Result of segmentation and selection by tak- 
ing into account a priori information on the center of the 
gravitational lens. The colors encode the objects eccen- 
tricity: [0.7,0.8] - green, [0.8,0.9] - yellow, [0.9,1.0] - 
red, i.e. the red objects are most likely arcs. The assumed 
center of the cluster is marked with a red cross. 




t • • 




Fig. 14. To illustrate the features of the new algorithm 
we compare some specific objects detected by both al- 
gorithms, the proposed algorithm (first row) and the seg- 
mentation by SExtractor (second row). The images of 
each column show the same part of the image. The pix- 
els detected by the according algorithm as belonging to 
the arc are plotted in red. In order to distinguish between 
close-by objects we use green color in addition. In col- 
umn 1 and 4 SExtractor has connected pixels to the 
arc which do actually come from various other sources. 
Therefore the arc may not satisfy the required elongation 
and hence may not be selected. Column 5 shows that the 
new algorithm has detected many more pixels of the faint 
structure than SExtractor. 
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Fig. 15. Result of segmentation and selection applied to 
the second test image. Parameter settings: Diffusion: filter 
parameter T - 30, edge sensitivity K - 0.0001, pre-filter 
parameter cr = 2 and pre-filter parameter for structure 
tensor - 30; Selection: intensity threshold for detec- 
tion Cmin = 0.1, minimal intensity difference from back- 
ground Cmin2 - 0.2, eccentricity threshold Cecc - 0.7 and 
thickness threshold CtMck = 9. Colors as in Fig. 1131 



Fig. 16. Result of segmentation and selection applied to 
the third test image. Parameter settings: Diffusion: filter 
parameter T - 20, edge sensitivity K - 0.0001, pre-filter 
parameter cr - 2 and pre-filter parameter for structure 
tensor yu = 15; Selection: intensity threshold for detec- 
tion Cmin = 0.3, minimal intensity difference from back- 
ground c,nin2 = 0.2, eccentricity threshold Cecc = 0.7 and 
thickness threshold Cthkk = 9. Colors as in Fig. 1131 
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3. Segmentation: the principal directions as determined in the 
diffusion process are exploited in the later segmentation 
process and the edge enhancing eases border allocation. 

4. Selection: after detection thin and elongated objects are se- 
lected. 

Comparing our algorithm with the software package 
"SExtractor" we find that both algorithms rely on multi-level 
strategies for feature extraction in astronomical data. The em- 
phasis in SExtractor Ues on the extraction of stars and galax- 
ies while our algorithm is designed to extract elongated arcs. 
Both multi-level strategies are implemented rather differently; 
the major differences are the filtering methods used and the de- 
termination of the segmentation threshold: SExtractor relies 
on Gaussian convolution filtering; the threshold for segmenta- 
tion depends on the estimated background. Our algorithm reUes 
on anisotropic diffusion filtering; the segmentation threshold is 
determined by border allocation using the same principal di- 
rections as in the diffusion process and avoiding the effect of 
merging close-by objects. 

The new algorithm is particularly well suited for the detec- 
tion of arcs in astronomical images. It can be both appUed to 
automated surveys as well as to individual clusters. 

The algorithm (implemented in C) will be pro- 
vided to the pubUc. Feel free to contact Frank Lenzen 
(Frank.Lenzen@uibk.ac.at). 
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